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PREDICTION  OF  POLYTOMOUS  EVENTS: 
MODEL  DESCRIPTION.  ALGORITHM  DEVELOPMENT 
AND  METHODOLOGICAL  ASPECTS.  WITH  AN  APPLICATION 

I.  O'MUIRCHEARTAIGH 
D.  P.  GAVER 

1 .   Introduction 

The  prediction  of  dichotomous  events  in  meteorology  (fog/no  fog, 
precipitation/no  precipitation)  has  been  widely  studied.   Such  predictions  are 
also  of  interest  in  reliability  and  survival  analysis,  and  in  manpower 
planning.   The  analysis  generally  involves  logistic  regression  or 
(equivalently)  linear  discriminant  functions.   Most  standard  statistical 
packages  (e.g.  BMDP,  SAS)  provide  the  facility  for  performing  this  analysis  in 
some  form.   Also,  the  book  by  Cox  (1970)  can  be  consulted. 

A  natural  extension  of  this  problem  (and  one  which  has  many  potential 
applications  in  meteorology  and  elsewhere)  is  the  situation  in  which  the 
predictand  is  polytomous,  i.e.  has  multiple  categories.   For  example,  it  might 
be  desired  to  predict  visibility  (good/marginal /bad)  or  precipitation 
( none/ rain/ snow ) .   Two  (methodologically)  distinct  cases  can  be  envisaged, 
viz. 

a.  when  the  predictand  involves  ordered  categories 

b.  when  the  categories  are  unordered. 

Typically,  the  former  case  is  the  more  common  in  meteorological  applications, 
and  this  is  the  problem  addressed  in  this  report. 

The  particular  application  analyzed  here  involves  583  records  of  time  to 
formation  of  tropical  storms,  and  associated  values  of  various  meteorological 


variables.   The  time  to  storm  formation  is  polytomized  and  recorded  as 

1:   storm  formed  within  24  hours 

2:   storm  formed  between  24  and  48  hours 

3:   storm  formed  between  48  and  72  hours 

4-*   storm  did  not  form 
The  five  meteorological  variables  recorded  were: 

X1 :   unconditional  probability  of  storm  formation  -  a  measure  of  the 
likelihood  of  storm  formation  in  the  particular  disturbance 
location  at  the  given  time  of  year. 

XQ:   large  scale  vorticity  (computed  over  5  Latitude  grid). 

X„:   divergence  (computed  over  5  Latitude  grid). 

X:   small  scale  vorticity  (computed  over  25  Latitude  grid) 

X'-      local  generation  of  vorticity  (product  of  X„  and  X4). 
Our  objective  was  to  determine  how  much  predictive  information  is  provided  by 
the  meteorological  variables  to  facilitate  prediction  of  imminence  of  tropical 
storms.   Essentially  the  problem  involves  regression  models  where  the 
dependent  variable  is  ordinal.   Much  attention  has  been  devoted  to  this 
general  problem  in  the  statistical  literature  of  recent  years  (McCullagh 
[1980],  McCullagh  and  Nelder  [1983],  Green  [1984],  Anderson  [1984]).   The 
central  concept  is  that  of  the  generalized  linear  model  (McCullagh  and  Nelder 
[1983]). 

In  section  2  we  describe  briefly  the  concept  of  generalized  models,  and 
in  more  detail,  the  particular  model  (an  extension  of  [dichotomous]  logistic 
regression)  utilized  for  our  data.   In  section  3  we  summarize  the  results  of 
an  ad  hoc  application  of  the  model  to  our  data,  present  the  relevant  parameter 
estimates,  and  evaluate  the  predictive  performance  of  our  model.   A  general 


discussion  of  our  results  is  presented  in  section  4,  together  with  suggestions 
for  future  work. 

2.   The  Model. 

2.1  General  formulation. 
Following  Green  [1984]  and  McCullagh  and  Nelder  [1983],  we  consider  a  log 
likelihood  L.  a  function  of  an  n-vector.  17,  of  predictors.   We  postulate,  in 
our  model,  that  the  predictor  17  is  functionally  dependent  on  the  p-vector  P   of 
parameters  of  interest.   For  our  particular  application,  17,  and  P   are 
specified  in  Section  2.2.   The  maximum  likelihood  estimation  of  P   involves 
essentially  solution  of  the  equations 


§=DTu  =  0  (1) 


where  u   ,  =  -r— 
nxl   017 


and  D    =|g 

nxp   6p 


using  the  notation  of  Green  [1984].   Again  following  Green,  the  standard 
Newton-Raphson  method  for  the  iterative  solution  of  (1)  involves  evaluating  u, 
D  and  the  second  derivatives  of  L  for  an  initial  value  of  /3,  and  then  solving 
the  linear  equations 


A)    (0*  -  P)   =   DTu  (2) 


for  an  updated  estimate  P   .      Green  shows  that  this  is  approximately  equivalent 


to  the  solution  of 


(DTAD)(P*-0)  =  DTu  (3) 


T 

where  A    =  E(-r—  (-r— )  ) 
nxn    votj  v6tj'  ' 


given  an  initial  estimate  of  P  (about  which  we  have  some  further  comments  in 
Section  2.3).  Equation  (3)  can  be  solved  directly  for  0  ,  or,  equivalently , 
if  a  weighted  least  squares  program  is  available.  p*  results  from  regression 
A  u  +  Dp   onto  the  columns  of  D  using  weight  matrix  A. 


2.2  Specific  Formulation. 
In  our  application  the  data  are  in  the  form  of  N  multinomial  samples  on 
the  same  set  of  k(=4)  response  categories  (e.g.  categories  1.  2,  3  and  4 
indicating  storm  imminence  as  described  in  Section  1).   The  data  may  be 
arranged  as  a  two-way  table  of  counts  y.  .,  i=l...N;  j  =  l,  ,,.4.   The  log- 
likelihood  L  is  then  given  by 


L  =  2  2  y   log  p  (4) 

i  j 


where  p.  .are  the  cell  probabilities,  and  2  p.  .  =  1 

1J  j-i   1J 


In  the  case  where  the  categories  1,2 k  are  ordered,  McCullagh  and 

Nelder  [1983]  and  Green  [1984]  both  suggest  the  model 


j         G.   2x.   /3 

p   =  +(-^— )  (5) 

J    r=l   ir         Ti 


i=l.  2 N;   j  =  1.2 k-1. 

where  n. .  represent  for  fixed  i,  the  cumulative  cell  probabilities,  the  matrix 
(x.  )  represents  covariate  information  and  ^   is  a  given  distribution  function. 
Motivation  for  their  model  is  provided  by  considering  the  response  variable  as 
an  arbitrary  grouping  of  an  unobservable  underlying  variable  on  a  continuous 
scale  with  "cutpoints"  6.,    .  ..8,  ...   In  some  applications  the  8  's  will  be 
unknown  and  will  need  to  be  estimated;  in  others,  such  as  the  present  one, 
they  will  be  known,  because  the  categorical  variable  (storm  imminence  coded 
1,2,3,4)  really  is  an  arbitrary  grouping  of  an  underlying  continuous  variable 
based  on  known  cutpoints  (in  this  case  that  variable  is  time  to  storm 
formation  with  cut-points  6^24,  8^=48  and  8„=72)  .   For  our  application,  we 
chose  \p   to  be  the  logistic  distribution  function,  viz.  , 


*(*)  =  -^  (6) 

1+e  X 


This  is  the  most  widely  used  model  in  applications  and  has  the  advantage  that 
a  simple  transformation  can  be  used  to  (a)  graphically  check  the  suitability 
of  the  model  and  (b)  provide  initial  values  for  the  iterative  estimation 
process. 


Our  model  therefore  is 


'/ 


ij  "  M 

e  -  i  x.  p 

J    m-1   lm  m 

,  m=l 

1  +  exp  -  

Ti 


(?) 


or.  in  terms  of  previous  notation 

T)  =  i7(P)  (8) 

where  /3  =  (/31  ,  .  .  ./3u,t1  ,  .  .  .tm)  .   Unless  we  impose  constraints,  identif  iabi  1  i  ty 
1      Ml      IN 

problems  can  arise.   One  common  expedient  (based  on  the  concept  of  an 
underlying  continuous  variable  used  for  the  classification)  is  to  allow  an 
intercept  and  scale  to  write 

M 

"ij  =  *< o^ )  <9> 

i  =  l....N;  j=l k-1. 


where  8's  are  known. 

A  reformulation,  more  convenient  for  actual  computation 


r,  =Wo  +  Ple  ♦  i,,^)  (io) 

m=l 


2.3  Specific  Methodology 
We  now  apply  the  general  methodology  of  Section  2.1  to  the  specific  model 


described  in  equation  (9).   This  involves  two  steps 

(a)  deriving  explicit  expressions  for  A,  D  and  u  described  in  that 
section,  and 

(b)  Finding  a  suitable  starting  value  P     for  the  iterative  reweighted 
least  squares  (IRLS)  procedure. 

We  describe  first  the  expressions  for  D,  A  and  u  for  the  special  case  of 
M=l  covariates  with  k=4  categories.   The  extension  to  other  cases  is 
straightforward.   If  we  let 

Vci  =  ^ir  T?i2f  ^n*  ^21 '  ^22 vm^ 


an 

where  n  =  Nx3.  then  D  =  ~tt  is  a  matrix  given  by 


D 


nx3 


FCBj.Xj) 
F^.x^ 
F(63.x1) 
F(61.x2) 


e1F(e1.x1) 
e2F(e2.Xl) 

9^(6^X2) 


XjFCBj.Xj) 
x1F(82,x1) 
x1F(G3.x1) 

x2F(erx2) 


F(93.xN)   e3F(e3,xN)   ^(83^) 


/ 


(11) 


where 


F(e..xj)=- 


exp  {Pj+  p^+^x.} 


Similarly,  u  ,    =  3—  is  given  by 
J        nxl   &q 


u 


11 


12 


^11   ^12  "  ^ll 


12 


13 


^"^ll   ^13  "  ^12 


13 


14 


^n-71^  l  ~  * 


13 


21 


22 


^21   ^22  "  ^21 


(12) 


N3 


N4 


^3  ^2   l     ^3 


and  A    =  F  ( —)    is  given  by 

nxn     v„ T'    &      J 

dim 


A  = 


'12 


T,11^12"T,11) 
1 
^12^11 


T)12-r?u) 
"(T,13"T?11) 


0 


0 


(T?12"T?11)(1713"T?12) 

1 


-d-n12) 


^lS"7^ 


Ko-lioK^ln) 


13  '12 


13' 


0 


-T?, 


22 


t701(t?00-t]01) 


'21 v  '22  '21 


etc 


(13) 


J 


0  0 

with  tridiagonal  3x3  matrices  similar  to  the  one  given  above  along  the  main 
diagonal,  and  zeroes  elsewhere. 


8 


Each  of  A,  D  and  u  depends  on  the  unknown  parameters  p*.   Given  an  initial 
value  for  p*  we  can  evaluate  A,  D  and  u  and  commence  the  iterative  process. 

The  initial  values  can  be  obtained  by  noting  that  if  we  apply  a  logit 
transformation  to  y.   in  equation  (10)  to  give 


tn  T^J  =  P°  +  piej  +  P2Xi  <14> 


then  these  logits  are  linear  in  the  parameters  g.   Initial  estimates  of  g  can 
be  obtained  by  constructing  empirical  logits. 


j         1 
^     yir>  +  2 
6.  ,  =  In  {±2— , }  (15) 


ij      l     j 


1 


n  -  (2  y.r  +  ^) 
r=l 


and  performing  an  unweighted  least  squares  analysis  for  the  model  of  equation 
(14)  using  these  empirical  logits  as  the  dependent  variable. 

Two  points  should  be  noted  in  relation  to  the  estimation  of  initial 
values  for  p*.   Firstly,  the  IRLS  procedure  is  quite  sensitive  to  the  choice  of 
initial  value,  (see  Green  (1984)),  and  a  poor  choice  can  lead  to 
non-convergence  of  the  algorithm.   Secondly,  to  obtain  initial  values  by  this 
procedure,  it  may  be  necessary  to  group  the  observations  into  categories  based 
on  values  of  the  independent  variables(s) .   If  the  data  are  not  so  grouped, 
(and  our  data  consists  essentially  of  multinomial  samples  of  size  1).  then  all 
y. .  will  be  either  0  or  1  in  our  case  this  will  lead  to  all  the  empirical 
logits  having  value  either  In  3  or  -In  3. 


Although,  it  may  be  necessary  to  group  the  data  to  obtain  starting  values 
for  P,    the  maximum  likelihood  estimation  of  P   may  be  carried  out  for  either 
the  grouped  or  the  ungrouped  data. 

2.4  Computational  procedure. 
Since  a  stepwise  program  was  not  available,  the  model  given  by  equation  (10) 
was  estimated  separately  for  each  covariate.   Using  the  deviance  (the 
likelihood  ratio  statistic  against  the  saturated  model)  as  a  measure  of 
goodness  of  fit.  the  best  single  explanatory  variable  was  X-.   Having  thus 
determined  the  optimal  single  variable  for  inclusion  in  the  model,  we  then 
proceeded  to  establish  which,  if  any.  variable  should  next  be  included  in  the 
model.   Due  to  the  non-availability  of  a  package  for  performing  this  analysis, 
the  computation  involved  was  cumbersome;  the  inclusion  of  an  additional 
variable  necessitated  the  re-programming  of  the  computations  leading  to  the 
matrices/vectors  A,  D  and  u.   It  was  determined  that  X,,  was  the  next  variable 
which  should  be  included  in  the  model.   A  third  step  of  the  stepwise  procedure 
was  also  carried  out.  but  no  additional  variable  warranted  inclusion  in  the 
model . 

Accordingly,  in  evaluating  the  predictive  performance  of  the  model,  and 
in  comparing  this  performance  with  those  of  discriminant  analysis  and  multiple 
regression,  we  used  only  the  explanatory  variables  X-   and  X~. 

3.   Evaluation  of  the  predictive  performance  of  the  model. 
3. 1   Introduction 
The  model  developed  in  this  paper  essentially  produces,  for  given  values 
of  the  covariates,  estimates  of  tj .  . ,  the  cumulative  category  probabilities, 


10 


and  from  these  estimates  of  p   ,  the  actual  category  probabilities.   Hence, 
the  model  provides  probability  forecasts  of  the  four  storm  imminence 
categories,  for  a  given  meteorological  situation  (as  represented  by  the  values 
of  the  meteorological  variables).   These  probabilistic  forecasts  can  be  given 
directly  as  such,  or  may  be  converted,  by  methods  described  in  Section  3.2. 
into  categorical  forecasts. 

The  problem  of  evaluating  statistical  forecasts  of  this  type  has  been, 
and  continues  to  be.  a  topic  of  major  interest  in  meteorology.   In  Section 
3.3,  we  consider  two  very  simple  methods  of  evaluating  such  forecasts;  these 
two  methods  do  not  necessarily  lead  to  the  same  conclusions  in  relation  to  the 
relative  performance  of  the  various  forecasts. 

3.2  Use  of  the  Model  for  Forecasting 
We  consider  two  possibilities: 

(a)  Given  the  estimated  probability  forecast,  a  categorical  forecast  can 
be  provided  by  forecasting  the  category  of  maximal  probability.   We  denote 
this  forecast  by  Fl . 

(b)  The  model  described  by  equation  (10)  (or,  more  intuitively  by 
equation  (9)),  suggests  an  underlying  continuous  variable,  say  Z,  with  the 
explanatory  variable  falling  into  categories  1.2,3  or  4  accordingly  as  Z  i  0^, 
Q     <  Z  <  e„.  92  <  Z  i   93.  Z  >  e3.  respectively.   Given  values  of  the 
covariates.  and  estimates  p*  of  p\  we  can  estimate  Z  by  (in  the  notation  of 
equation  (10)) 
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.   M   „ 

Ho        .    nTm+2 
Z  = E^ (17) 

Pi 


and  then  provide  a  categorical  prediction  that  the  storm  imminence  category  is 
1,  2,  3  or  4  according  as 


z  *  8j  9j<  z  *  e2.  e2  <  z  s  e3.  z  >  e3 


This  forecast  is  denoted  by  F2. 

3.3  Evaluation  of  the  forecasts 

Since  the  two  forecasts  which  we  are  considering  here  are  categorical 
forecasts,  one  plausible  criterion  for  evaluating  these  forecasts  would  be  the 
number  of  correct  forecasts  achieved.   Different  forecasts  can  be  readily 
compared  using  this  measure.   Since  the  categories  being  forecast  are  ordered, 
an  incorrect  forecast  which  is  within  one  category  of  being  correct  is 
presumably  preferable  to  one  which  "misses"  by  two  or  more  categories.   Hence 
an  alternative  measure  of  performance  would  be  the  number  of  forecasts  which 
are  within  one  category  of  being  correct.   We  use  both  of  the  above  measures 
in  the  paper  to  compare  forecasts.   As  we  will  show,  the  different  criteria 
can.  in  some  cases,  lead  to  a  different  ranking  of  forecasts. 

In  estimating  the  predictive  performance  of  the  model  using  the  above 
measures,  we  omitted  each  data  point  in  turn,  estimated  the  model  parameters 
from  all  the  remaining  data,  and  then  used  the  forecast  procedures  Fl  and  F2 
to  predict  the  category  of  the  omitted  data  point.   For  comparative  purposes, 
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we  also  used  the  standard  techniques  of  discriminant  analysis  and  multiple 
regression  (using  in  the  latter  case  as  dependent  variable  the  coded  imminence 
of  tropical  storm  variable,  which  takes  values  1.  2.  3  and  4. 

4.   Discussion 

The  results  of  the  cross-validation  procedure,  described  in  Section  3.3. 
for  forecast  methods  Fl  and  F2.  and  for  discriminant  analysis  and  multiple 
regression,  are  given  in  Tables  1.2,3  and  4  respectively.   An  overall  summary 
of  the  relative  performance  of  the  various  methods  is  presented  in  Table  5. 
We  would  emphasise  that  any  conclusions  drawn  here  in  relation  to  the  efficacy 
of  the  various  procedures  are  valid  only  in  relation  to  the  present 
application.   Broader  statements  about  the  general  performance  of  these 
methods  would  require  extensive  further  analysis. 

It  is  clear  from  Table  5  that  no  single  technique  is  clearly  superior. 
Using  the  criterion  of  maximising  the  number  of  correct  forecasts,  the 
generalized  linear  model  applied  in  this  paper,  with  forecasting  strategy  F2, 
is  the  best  among  those  considered.   However,  if  maximising  the  numbers  of 
forecasts  correct  within  one  category  is  chosen  as  the  comparative  criterion, 
then  multiple  regression  emerges  as  the  optimal  methodology. 

From  Table  5  it  is  clear  that  the  performance  of  discriminant  analysis 
is.  in  this  application  at  least,  somewhat  inferior  to  that  of  the  other 
techniques.   A  comparison  of  Tables  1  and  2  (which  use  the  same  model  but 
different  forecasting  strategies  based  on  the  model)  reveals  some  interesting 
facets  about  each  of  these  strategies.   Procedure  Fl  always  forecasts  either 
category  1  or  category  4,  and  produces  the  highest  overall  percentage  of 
correct  forecasts.   Procedure  F2  is  less  extreme,  and  "smears"  the  category  1 
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forecasts  over  categories  1  and  2,  and  the  category  4  forecasts  over 
categories  3  and  4;  this  has  the  effect  of  reducing  the  percentage  of  correct 
forecasts  but  increasing  the  percentage  of  forecasts  correct  to  within  one 
category.   Multiple  regression  "smears"  the  forecasts  still  further,  with  a 
resultant  decrease  in  the  percentage  of  correct  forecasts  and  increase  in  the 
percentage  correct  to  within  one  category.   Multiple  regression  is,  in  fact, 
the  method  which  produces  the  highest  overall  percentage  of  forecasts  correct 
to  within  one  category. 

It  is  clear  that  the  choice  of  strategy  (forecasting  procedure)  among 
those  considered  and  described  here  will  be  greatly  influenced  by  the  relative 
importance/seriousness  of  the  various  correct/incorrect  forecasts.   This 
strongly  suggests  that  a  decision  theoretic  approach  might  be  considered, 
although,  in  practice,  the  specification  of  a  loss  function  may  be  difficult 
and  may  unduly  influence  the  choice  of  strategy. 
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Table  1 
Forecast  Procedure:   Fl 


Number  of  Forecasts 

True 
State 

Forecast  State 

1 

2 

3 

4 

Total 

1 

27 

0 

0 

34 

61 

2 

10 

0 

0 

32 

42 

3 

4 

0 

0 

20 

24 

4 

11 

0 

0 

445 

456 

total 

52 

0 

0 

531 

583 

Table  2 
Forecast  Procedure:   F2 


Number  of  Forecasts 

True 
State 

Forecast  State 

1 

2 

3 

4 

Total 

1 

18 

9 

7 

27 

61 

2 

8 

2 

5 

27 

42 

3 

1 

3 

3 

17 

24 

4 

4 

7 

17 

428 

456 

total 

31 

21 

32 

499 

583 
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Table  3 
Forecast  Procedure:   Discriminant  Analysis 


Number  of  Forecasts 

True 
State 

Forecast  State 

1 

2 

3 

4 

Total 

1 

27 

12 

16 

6 

61 

2 

14 

5 

9 

14 

42 

3 

3 

3 

14 

4 

24 

4 

14 

26 

62 

354 

456 

total 

58 

46 

101 

378 

583 

Table  4 
Forecast  Procedure:   Multiple  Regression 


Number  of  Forecasts 

True 
State 

Forecast  State 

1 

2 

3 

4 

Total 

1 

4 

21 

31 

5 

61 

2 

0 

11 

21 

10 

42 

3 

0 

2 

18 

4 

24 

4 

0 

11 

110 

335 

456 

total 

4 

405 

180 

354 

583 
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Table  5 


Forecast 
Procedure 

%  Correct 
Forecasts 

%  Forecasts 
Correct  with 
one  Category 

Fl 

81 

86 

F2 

77 

88 

Discriminant 
Analysis 

69 

86 

Multiple 
Regression 

63 

90 
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